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Abstract 

The role of a matrix response to a fluid insertion is analyzed in terms 
of a perturbation theory and Monte Carlo simulations applied to a hard 
sphere fluid in a slit of fluctuating density-dependent width. It is demon- 
strated that a coupling of the fluid-slit repulsion, spatial confinement and 
the matrix dilatation acts as an effective fluid-fluid attraction, inducing a 
pseudo-critical state with divergent linear compressibility and non-critical 
density fluctuations. An appropriate combination of the dilatation rate, 
fluid density and the slit size leads to the fluid states with negative linear 
compressibility. It is shown that the switching from positive to negative 
compressibility is accompanied by an abrupt change in the packing mecha- 
nism. 

I. INTRODUCTION 

Compounds with a negative linear (or surface) compressibility have recently attracted 
an interest 1 3 because of their specific properties (stretch-induced densification and aux- 
etic behavior), which might have promising applications. In the case of pure materials only 
"rare" crystal phases exhibit these effects, while for composite structures 4 it seems to be 
rather common. Recent experimental studies on insertion into organic 5 and non-organic 6 
matrices reveal the negative compressibility effects due to the host-guest coupling. Con- 
ceptually similar escape transition 7 occurs when a polymer chain is compressed between 



two pistons. Charge separation in confined fluids has also been accompanied 8 with the 
negative compressibility. These examples suggest that the insertion systems with confine- 
ment and swelling (dilatation) should exhibit this generic effect at appropriate conditions. 
In particular, two-particle systems confined to two-dimensional finite-size boxes of differ- 
ent geometry have been studied 9-11 in this context. It has been found that the isotherms 
exhibit a van der Waals-type (vdW) instability (loop) as the box size passes through 
a critical value. The instability has been associated with a prototype of a liquid-gas 
or a liquid-solid transition in many-body systems. A quite similar instability appears in 
two-dimensional granular media 12 , whose effective temperature decreases with the density, 
leading to a phase separation. Similar features have been detected 13 in hard discs confined 
to a narrow channel. Nevertheless, the non-monotonic behavior has been attributed to a 
sharp change in the accessible phase space, without any connection to collective effects 
typical for the conventional phase transitions. It is well-known, however, that a liquid-gas 
transition in finite systems 14 manifests itself through the negative compressibility states 
due to the surface effects. 

Thus, it is quite interesting to find out if such a loop could exist in three-dimensional 
many-body systems and what is the physics behind. In particular, we focus on a proto- 
type of an insertion system. One of the key features of such systems is a host dilatation 
(or contraction) upon a guest accommodation. This is evident from experiments with 
high-porosity materials, like aerogels 15 ' 16 . Upon adsorption such matrices change in vol- 
ume and their pore size distribution depends on the adsorbate pressure (or density). This 
effect is not exclusive to relatively soft gel-like matrices, it is quite common for carbon 
nanotubes 17 and various intercalation compounds 18-20 . In anisotropic cases (e.g. layered 
compounds) one deals with a competition of two effects. An increase of the lateral dimen- 
sion (stretching) at constant number of particles tends to decrease the guest density. This 
usually induces a transversal shrink, leading to a densification in this direction. Since we 
have a composite host-guest system, the shrink does not obey the linear elasticity rules 
and, depending on the shrink intensity, one might expect the negative compressibility. 
In particular, it has been shown 20 that a nonlinear increase of the host size with the 



guest concentration may induce a liquid-gas coexistence for the guest fluid even if the 
bulk liquid phase does not exist (e.g. a hard-sphere fluid). This means that the negative 
compressibility states were indeed present, but they were unstable under the conditions 
imposed. 

Therefore, the objective of this paper is to find the conditions at which a coupling of the 
host dilatation and the guest confinement can stabilize the negative compressibility states. 
For this purpose we start with a quite simple model - the hard spheres, adsorbed in a softly- 
repulsing planar slit with a density-dependent width. The model contains all the relevant 
features: confinement, dilatation and spatial anisotropy. The latter is important because 
it offers a possibility of observing negative linear and positive tangential compressibilities, 
preserving the overall thermodynamic stability. On the other hand, in the absence of the 
dilatation (fixed slit width), the model is phenomenologically simple 21,22 . In particular, 
there is no a liquid- vapor transition (both in the bulk and in confined geometries). Such 
that the new aspects are not masked by the internal complexity. The system is analyzed 
by means of Monte Carlo (MC) simulation and a first-order perturbation theory 23-25 . 
Qualitative reliability of this approach has been tested in application to the liquid-vapor 
coexistence 24 and phase separation 25 in confined geometries. 



Consider N hard spheres of diameter a = 1 in a planar slit with the surface area S. 
The slit itself is a part of a hosting system, whose coupling to the guest species changes 
the slit geometry. Therefore, the slit width h is not fixed, it fluctuates according to the 
fluid density (see below). For any h the total Hamiltonian is 



II. MODEL 



H = H ff + H f 



(1) 



where Hff is the hard sphere Hamiltonian, Hf w is the slit potential 



1 



1/2 < Zi < h - 1/2. 



(2) 



i=l l^i 



(h - Zi ) k 



We do not take into account a short-ranged fluid-wall attraction, responsible for the 
surface adsorption or layering effects. The inverse-power shape for Hf w is chosen as a 
generic form of a soft repulsion, with k controlling the softness. For technical purposes we 
are working with k — 3. Moreover, our results are qualitatively insensitive to a particular 
choice of k. 

It should be noted that we do not discuss an adsorption mechanism or the equilibrium 
between the pore and the bulk fluids. In our case N is fixed and we focus on the pressure 
variation due to the changes in the slit geometry. 



III. PERTURBATION THEORY 

In practice the evolution of the matrix morphology is much slower than the fluid 
equilibration process. Then for a given width h we can calculate the fluid thermodynamics 
conditional to h. 



A. Conditional equation of state and insertion isotherm 

At any pore width h the free energy can be represented as 

PF(h) = /3F {h) - ln(e-^) (3) 

where F (h) is the free energy of a reference system, and (...) is the average over the 
reference state. Taking a spatially confined hard sphere system (the one with A = 0) as 
a reference, we consider a first order perturbation 23-25 for the conditional free energy 

(3F(h) = (3F (h) + pSp(h) f [] dziH fw (4) 

where (3 = l/(kT) and p(h) is the pore density in the "slab" approximation 

N 

p{h) = W^Ty (5) 

ignoring a non-monotonic behavior of the density profile with increasing pore density. 
This is reasonable for wide pores and low fluid densities. The reference part is estimated 



in the excluded volume approximation, while the perturbation contribution is N^(h), 
such that the total conditional free energy is 

1 - bp(h)~ 



pF(h) = -Nln 



N + N[3^(h) (6) 



A* T p(h) 

where b = 2na 3 /3 is the excluded volume factor, A T is the thermal de Broglie length and 
Tangential Pt(h) and normal P n (h) pressures can be found as 



This leads to 



fl p/n JL. BP(h ) PW i l6A * N 2h+1 (a) 

where A* = (3 A. Equation (5) allows us to eliminate the surface density N/ S in the favor 
of p{h) and we obtain the following equation of state 

Therefore, for a fixed h, our result is clear and simple. The tangential pressure has the 
bulk form, with the bulk density being replaced by the pore density. As a consequence of 
the first-order perturbative approach, Pt(h) does not depend on the slit-fluid interaction. 
As we will see later, the simulation results demonstrate that this dependence is indeed 
quite weak. If necessary this effect can be reproduced theoretically taking into account 
the second-order perturbation term. The normal pressure increases due to the repulsion 
A* and decreases with increasing pore width h, such that P n (h) = Pt(h) as h — > oo. 

Recall that we are dealing with a fixed N. The insertion process can be considered 
in the same framework. Just instead of the pressure components one would calculate the 
system response to increasing N- the insertion isotherm 

A 3 T p(h) 



V uiy / p,s,h 



l-bp(h) 

which relates the fluid density and the chemical potential p{h). 



B. Dilatation effect 



As is discussed above, the pore dilatation can be taken into account, assuming that 
the width h is density-dependent. Real insertion materials are usually rather complex 
(multicomponent and heterogeneous). For that reason one usually deals with a distribu- 
tion of pore sizes or with an average size. In this context we assume that h is known only 
statistically and the probability distribution f(h\p) is conditional 26 ' 27 to the guest density 
p, which should be selfconsistently found from 

p = J dhf(h\p)p(h). (12) 

Then, taking our results for P n (h) and Pt(h), we can focus on the equation of state 
averaged over the width fluctuations. 

Pi = j dhf(h\p)Pi(h); i = n,t. (13) 

This, however, requires a knowledge on the distribution f{h\p). Even without resorting 
to a concrete form for f(h\p), it is clear that the matrix reaction can be manifested as 
a change in the distribution width or/and the mean value. One of the simplest forms 
reflecting at least one of these features is a 5-like distribution, ignoring a non-zero width. 

f(h\p)=5[h-h(p)} (14) 

where S(x) is the Dirac 5-function and the mean pore width h(p) is density-dependent. 

h(p) = h (l + tanh[A(p - p )] + tanh[Ap ]) (15) 

This form mimics a non-Vegard behavior, typical for layered intercalation compounds 20 . 
At low densities (p << p ) the dilatation is weak. The most intensive response is at 
p ~ po, and then the pore reaches a saturation, corresponding to its mechanical stability 
limit. Here A is the matrix response constant or dilatation rate, controlling the slope 
near p m p . 

From eq. (12) the average density is found to be 



Changing the surface density N/ S we vary the average pore density p. This allows us to 
eliminate N/S in the favor of p in all thermodynamic functions. Combining eqs (10) and 
(13) we obtain 

^ = 1^ + 1M P (2%) - 1)3 ' (17) 

It is convenient to introduce the compressibility function 

i dhjp) _ i dh{ P ) dp 
Xn h(p) - 1 dP n h( P ) - 1 dp dP n 1 ' 

which vanishes in a non-swelling limit (h(p) = h ). It is clearly seen that Xn could be 
singular either if the swelling is discontinuous, or if P n exhibits an inflection point. In 
what follows we focus on the second option. 

In the case of a wide and weakly reacting pore we expand in terms of 1/ho and A, 
obtaining a generic van der Waals form 

1 — bp ho h 

Therefore, we have an interplay of several effects - the packing (first term), the fluid- 
matrix interaction (linear in density), and the matrix reaction (quadratic term). It is 
seen that a coupling of the fluid-slit repulsion (A), spatial confinement (ho) and the 
matrix response (A) acts as an effective infinite-range fluid-fluid attraction (but just in 
one direction). Introducing a dimensionless temperature T* = 1/(4A*), and solving 

dP n d 2 P n , s 
we find the "critical" parameters 

+ r l (2A-6)3 

Pc "3 Ab ' c 27 bhoA* 1 } 

at which the normal compressibility Xn diverges. It is seen that a physically meaningful 
(with T* > 0) pseudo-criticality appears only at A > 6/2. In other words, the matrix 



reaction A should dominate the packing effects. In addition, T* decreases with increasing 
pore width h . 

It is well known that for bulk systems the vdW loop appearing at T* < T* is un- 
physical and one usually invokes the Maxwell construction, determining the liquid-vapor 
coexistence. As we will see below, in our case the loop is a physically justified effect. It 
appears as a competition between the packing, which tends to increase the pressure with 
increasing p and the slit dilatation, decreasing P n with increasing p. As a result the states 
with negative linear compressibility are stabilized. Note that this behavior is not sensitive 
to the particular form (15) of h(p). Moreover, the dilatation rate A is essential, while the 
dilatation magnitude h(p) — h plays only a marginal role. As it should be, this effect 
disappears in the bulk limit ho — > oo or in the case of insensitive matrices A — > 0. 

In order to be more accurate in comparison to the MC data the reference part is 
replaced by the Carnahan-Starling form 

+ (2%) + l)(%)-l) 
P p n = P {1 _ v]3 + p (2%) _ 1)3 , V = -P/6. (22) 

This allows us to avoid the unphysical behavior at high densities. 



IV. SIMULATION 

In order to verify the existence of the loop predicted by the perturbation theory as 
well as to get more insight into the phenomenon the MC simulation of the model has been 
carried out. We applied the canonical NVT MC simulations of a confined hard sphere 
fluid to calculate the density profiles and pressure components. The simulation cell was 
parallelepiped in shape, with parallel walls at surface separation h, and surface area 
S = L x x L y . The periodic boundary conditions were applied to the X and Y directions 
of the simulation box; the box length in the Z direction is fixed by the pore width. For 
a given pore width the adsorbed fluid density is chosen according to the Eq.(5). There 
was constant number of particles, N = 1000, and a desired fluid density has been get 
by adjusting the value of area, S. We repeated our simulation runs with bigger number 



of particles, N = 1500 and 3000, but no significant differences were found. The density 
profiles, p(z), were calculated in the usual way by counting the number of particles, iVj 
in the slabs of thickness d z = 0.05 parallel to the plane XY by using p(z) = Ni/v, where 
v is the slab volume v = d z L y . The definition of Irving and Kirkwood 28 has been 

used to calculate the components of the pressure. The normal component of the pressure 
for the fluid-fluid interaction is 

PPM = M - 1 (e E ^ ie(^)9(^)) (23) 

° \ i j>i ar Zt i z% i I 

where p(z) is the total density profile and (...) denotes the average over the MC configu- 
rations. The tangential component, Pt(z), has been calculated by substitution of zf^ by 
0.5(xjj + yfj) in equation (23). For hard spheres the derivative of the potential is 

= -5(r - 1) (24) 

dr 

The Dirac 5 function in our simulation was approximated as 

Q( r - 1) - Q(r - 1 - A) 
5(r - 1) = ^ L, (25) 

as A — > 0. In Eq. (25) A is used to define the region where two particles collides. A 
collision between particles occurs if 1 < r < 1 + A. Therefore, for the particular case 
of hard spheres the averaged normal component of the pressure tensor for fluid-fluid 
interactions can be reduced to 

v \ i j>i il I' U'l I 

The wall-fluid interaction is a function of z and affects only P n while the tangential 
component remains the same. The normal pressure in this case is 

PP?(z) = /3P?\z) + (3Pr(z) (27) 

where f3P^ /1 {z) and f3P^ 2 (z) are the fluid-wall contributions from surfaces located at 
zwi = 0.5 and zw2 = h — 0.5, respectively. These contributions are defined as 

PP?\z) = "f (E ^£^1* - *«H e (* - *) e (* - z ^)) ( 28 ) 




dH) w (z) 



dz 



ZW2 ~ Zi\Q(z W2 ~ Z)Q(Z - Zi) 



(29) 



These definitions of the pressure tensor treat the wall-particle interaction as a contri- 
bution to the intermolecular forces in a system consisting of the fluid and solid, rather 
than as an external field acting on the fluid. Thus the normal component of the pressure 
tensor must be independent of z as a condition of mechanical equilibrium and for a sys- 
tem containing bulk fluid must be equal to the pressure at the bulk density. Tangential 
components should approach the bulk pressure for sufficiently large h. 

In this work we have taken parameter A equal to 0.001, 0.002, and 0.003 to make the 
extrapolation. The average normal component, includes the fluid-fluid and fluid- wall con- 
tributions given by equations (26) and (27). Each simulation runs 5 x 10 5 MC cycles, with 
the first half for the system to reach equilibrium whereas the second half for evaluating 
the ensemble averages. 



In agreement with our theoretical prediction, the simulation results confirm that the 
negative compressibility states (the loop) appear only if the slit reaction A reaches some 
threshold value A* which involves a combination of h , p and A. The variation of the 
average density p due to the dilatation should dominate its variation, induced by the 
changes in the surface density N/S. 

Pressure as a function of the average pore density is plotted in Figure 1. It is seen 
that the normal component P n develops the loop with increasing pore-fluid repulsion 
A. The simulation results confirm that this feature is not an artefact of our theoretical 
approximations. As expected 25 , the perturbation theory systematically underestimates 
the magnitude of the fluid- wall repulsion, A, (overestimates the temperature T*) at which 
this effect takes place. The tangential component P t is much less sensitive to the slit 
dilatation. This is coherent with our theoretical estimation. Interestingly, that P t can 
be reasonably fitted by the expression for P n taken at much lower A (see the insets). 



V. RESULTS 



Therefore, we expect that at sufficiently strong repulsion (low temperatures) P t could 
also be non-monotonic. This would correspond to a fluid-solid transformation, which is 
not considered here. The loop becomes more pronounced with increasing dilatation rate 
A. Since only one pressure component exhibits this behavior, there is no rational basis 
to suspect a vdW instability (of liquid- vapor type). Moreover, we are dealing with rather 
wide pores (ho = 10) in order to avoid the narrow channel effects 13 . This makes us to 
search for an alternative explanation. 

With this purpose we have analyzed the fluid structure at different densities. The 
density profiles are presented in Figure 2. It is obvious that the system does not exhibit 
strong density oscillations. It is worth noting, that our perturbation theory results are 
obtained in the 'slab' approximation, eq.(5), which does not take them into account. On 
the other hand, analysis of the density profiles indicates a sudden change of the packing 
regime in the negative compressibility region (in between the points A and B, marked in 
the Figure 1 (b)). The fluid film becomes more dilute with decreasing (e. g. due to a 
lateral stretching) density up to the point B with p = 0.32. This process goes uniformly 
in the middle of the film and at the periphery. Passing from p = 0.32 to p = 0.29 does 
not change the middle density, while the rarefaction takes place only near the slit walls. 
Upon reaching p = 0.25 the film suddenly densities in the middle. This corresponds to 
the inflection point in Figure 1 (b). Up to the point A with p = 0.2 the film dilutes only 
at its periphery. Then we return to the usual uniform dilution with further decrease in 
the density. Therefore, there is a clear correlation between the negative compressibility 
states and the packing mechanism. 

In order to study the redistribution of the fluid density we have calculated the average 
density T(z) in slabs of different thickness z 



The obtained results are presented in Figure 3. As seen, F(z) exhibits the same loop as 
the normal pressure does. Going towards the middle of the pore makes this effect more 
pronounced. Such similarity permits us to conclude that the fluid-wall soft repulsion 




(30) 



(coupled to the dilatation A) is a main reason of the fluid reordering inside the pore, and 
as a consequence the states with negative compressibility can occur. This is clearly seen 
from eq. (22), that gives a monotonic isotherm as A — > or A — > 0. 

It is interesting to mention that a similar trend has been described recently in the 
studies of mineral clays swelling 30,31 . Namely, Smith et al. 31 reported the simulation 
results of hydrated Na-smectites with variable layer charge. They found that the tendency 
to swell increases with increasing layer charge (increasing repulsion), which is consistent 
with our conclusions. 

VI. CONCLUSION 

We have found that a coupling of the fluid-slit repulsion, spatial confinement and the 
slit dilatation acts as an anisotropic fluid-fluid attraction, inducing negative linear com- 
pressibility states. These states are shown to be related to an abrupt change in the packing 
regime, including a local densification in the middle with decreasing surface density. This 
resembles the stretch-densification effects 1 ' 2 in negative compressibility materials. In the 
context of our model the mechanism is the following. Stretching in the lateral direction 
decreases the surface N/S and the pore p densities. This leads to a transversal shrink 
h(p) at the rate A. This stabilizes the negative transversal compressibility when the rate 
passes some threshould value. This mechanism is quite different from that explored in the 
low-dimensional systems 9-11 ' 13 , where the loops appeared essentially due to the small-size 
effects restricting the phase space accessibility when the box length became comparable 
with the particle size. As a result, the vdW feature is found 11 to be very sensitive to the 
box geometry (rectangular or spherical). In our case we deal with a three-dimensional 
system where the particles can exchange their positions almost freely (10 < h(p)/a < 30), 
although the phase space is somewhat restricted by the slit-fluid repulsion. Nevertheless, 
our model shares some of the low-dimensional features 9-11 ' 13 and we recover the usual 
monotonic behavior in the bulk limit ho — > oo. 

Since the tangential pressure component does not exhibit this effect, the system does 



not undergo a phase transition (at least in its traditional sense). This is confirmed by the 
absence of strong density fluctuations (in both directions). On the other hand, the loop 
appears in the direction, along which our system is finite (finite h) and, therefore, the 
negative compressibility can be considered as a finite-size effect 14 , vanishing in the bulk 
limit. Nevertheless, an inclusion of an attractive fluid- fluid interaction would result 23,24 
in the pore condensation. In this respect it would be quite interesting to analyze how 
the above stretch-densification coexists with the true critical behavior. The role of more 
isotropic geometries (e.g. spherical pores) as well as that of a non-zero distribution width 
(f(h\p)) and the potential softness k could also be discussed. Another interesting point is 
to describe the sorption behavior, that is, changing p by varying N at fixed S. In this way 
we can mimic a "dosen" adsorption 29 or the equilibrium with an infinite bulk reservoir. 
We plan to analyze these issues in a future work. 
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FIGURES 

FIG. 1. Normal and tangential pressures (the insets) in the case of a weakly (a) and strongly 
(b) repulsive slit. The other parameters are ho = 10, A = 15, po = 0.3 

FIG. 2. MC density profiles at different densities (indicated). Other parameters as in the 
previous figure. 

FIG. 3. Average fluid density in slabs of different thickness z (lines) and the normal pressure 
component (line with symbols). 
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